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Discretization of the viscous terms in current finite-volume unstructured-grid schemes are compared us- 
ing node-centered and cell-centered approaches in two dimensions. Accuracy and efficiency are studied for 
six nominally second-order accurate schemes: a node-centered scheme, cell-centered node-averaging schemes 
with and without clipping, and cell-centered schemes with unweighted, weighted, and approximately mapped 
least-square face gradient reconstruction. The grids considered range from structured (regular) grids to irreg- 
ular grids composed of arbitrary mixtures of triangles and quadrilaterals, including random perturbations of 
the grid points to bring out the worst possible behavior of the solution. Two classes of tests are considered. The 
first class of tests involves smooth manufactured solutions on both isotropic and highly anisotropic grids with 
discontinuous metrics, typical of those encountered in grid adaptation. The second class concerns solutions 
and grids varying strongly anisotropically over a curved body, typical of those encountered in high-Reynolds 
number turbulent flow simulations. Results from the first class indicate the face least-square methods, the 
node-averaging method without clipping, and the node-centered method demonstrate second-order conver- 
gence of discretization errors with very similar accuracies per degree of freedom. The second class of tests 
are more discriminating. The node-centered scheme is always second order with an accuracy and complexity 
in linearization comparable to the best of the cell-centered schemes. In comparison, the cell-centered node- 
averaging schemes are less accurate, have a higher complexity in linearization, and can fail to converge to the 
exact solution when clipping of the node-averaged values is used. The cell-centered schemes using least-square 
face gradient reconstruction have more compact stencils with a complexity similar to the complexity of the 
node-centered scheme. For simulations on highly anisotropic curved grids, the least-square methods have to 
be amended either by introducing a local mapping of the surface anisotropy or modifying the scheme stencil 
to reflect the direction of strong coupling. 


I. Introduction 

Both node-centered and cell-centered finite-volume discretizations are widely used for complex three-dimensional 
turbulent simulations in aerospace applications. The relative advantages of the two approaches have been extensively 
studied in the search for methods that are accurate, efficient, and robust over the broadest possible range of grid and 
solution parameters. The topic was discussed in a panel session at the 2007 AIAA CFD conference, but a consensus 
did not emerge. One of the difficulties in assessing the two approaches is that comparative calculations were not 
completed in a controlled environment, i.e., computations were made with different codes and different degrees of 
freedom and the exact solutions were not known. 

In this paper, we provide a controlled environment for comparing a subset of the discretization elements needed 
in turbulent simulations, namely that of the viscous discretization. In particular, we consider only Poisson’s equation 
as a model of viscous discretization. We use the method of manufactured solution, so that the exact solution is 
known, and conduct theoretical and computational studies of both accuracy and efficiency for a range of grids. There 
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are six main schemes considered: a node-centered scheme, cell-centered node-averaging schemes with and without 
clipping, and cell-centered schemes with unweighted, weighted, and approximately mapped least-square face gradient 
reconstruction. Each of the schemes considered are nominally second-order accurate. 

The two-dimensional (2D) grids considered range from structured (regular) grids to irregular grids composed of 
arbitrary mixtures of triangles and quadrilaterals. Highly irregular grids are deliberately constructed through random 
perturbations of structured grids to bring out the worst possible behavior of the solution. Two classes of tests are 
considered. The first class of tests involves smooth manufactured solutions on both isotropic and highly anisotropic 
grids with discontinuous metrics, typical of those encountered in grid adaptation. The second class of tests concerns 
solutions and grids varying strongly anisotropic ally over a curved body, typical of those encountered in high-Reynolds 
number turbulent flow simulations. The properties to be compared are (1) computational complexity (operation count) 
at equivalent discretization accuracy and (2) efficiency and robustness of iterative solvers. 

This paper is the first part from an envisioned series of three papers comparing node centered and cell centered dis- 
cretizations. The second part will address inviscid fluxes, and the third part will discuss complete Reynolds-Averaged 
Navier-Stokes formulations. 


II. Grid specifi cation 

This paper studies finite-volume discretization (FVD) schemes for viscous fluxes on grids that are loosely de- 
fined as irregular. There is no commonly accepted definition for irregular grids, so, for clarity, we specify the grid 
terminology used in this paper. 

A grid is classified as periodic if it has (1) a periodic node connectivity pattern (i.e., the number of edges per 
node changes periodically) and (2) a periodic cell distribution (i.e., the grid is composed of periodically repeated 
combinations of cells). Thus, periodic grids can be analyzed by Fourier analysis. Grids that are derived from periodic 
grids by a smooth mapping are called regular grids. Examples of regular grids include, but are not limited to, grids 
derived from Cartesian ones: triangular grids obtained by diagonal splitting with a periodic pattern, smoothly stretched 
grids, skewed grids, smooth curvilinear grids, etc. Grids that cannot be smoothly mapped to a periodic grid are called 
irregular grids. Grids with varying local topology are called unstructured, e.g., grids with the number of edges 
changing from node to node with no pattern. 

The regular and irregular grids considered in this paper are derived from an underlying (possibly mapped) Cartesian 
grid with meshsizes h x and h y and the aspect ratio A = £“■; both meshsizes of the underlying grid are assumed to be 
small, h y <C 1, h x <C 1. Irregularities are introduced locally and do not affect grid topology and metrics outside of a 
few neighboring cells. A local grid perturbation is called random if it is independent of local perturbations introduced 
beyond some immediate neighborhood. For computational grids generated for the reported studies, grid irregularities 
are introduced in two ways (both local and random): (1) the quadrilateral cells of the underlying grid are randomly 
split (or not split) into triangles; (2) the grid nodes are perturbed from their original positions by random shifts; the 
shifts are fractions of a local meshsize. 

Six grid types are considered: (I) regular quadrilateral (i.e., mapped Cartesian) grids; (II) regular triangular grids 
derived from the regular quadrilateral grids by the same diagonal splitting of each quadrangle; (III) random triangular 
grids, in which regular quadrangles are split by randomly chosen diagonals, each diagonal orientation occurring with 
probability of half; (IV) perturbed triangular grids, which are random triangular grids, with grid nodes perturbed from 
their initial positions by random shifts; (V) perturbed quadrilateral grids, which are quadrilateral grids with randomly 
perturbed grid nodes; and (VI) perturbed mixed-element grids, in which perturbed quadrangles are randomly split or 
not split by randomly chosen diagonals. Grids of types (III), (IV), and (VI) are irregular because there is no periodic 
connectivity pattern. Grids of types (IV)-(VI) are irregular because there is no periodic pattern for cell distribution. 
The representative grids are shown in Figure 1 . 

Our main interest is in studying accuracy and stability properties of FVD schemes on general irregular (mostly 
unstructured) grids with a minimum set of constraints. In particular, we do not require any grid smoothness, neither 
on individual grids nor in the limit of grid refinement. The only major requirement for a sequence of refined grids 
is to satisfy the consistent refinement property. The property requires the maximum distance across the grid cells to 
decrease consistently with increase of the total number of grid points, N. In particular, the maximum distance should 
tend to zero as N -1 / 2 in 2D computations. For three-dimensional (3D) unstructured grids, the consistent refinement 
property is studied in Ref. 21 

The locations of discrete solutions are called data points. For consistency with the 3D terminology, the 2D cell 
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(a) Type (I): Regular quadrilateral grid (b) Type (II): Regular triangular grid (c) Type (III): Random triangular grid 



(d) Type (IV): Perturbed triangular grid (e) Type (V): Perturbed quadrilateral grid (f) Type (VI): Perturbed mixed grid 


Figure 1. Typical regular and irregular grids. 


boundaries are called faces, and the term “edge” refers to a line, possibly virtual, connecting the neighboring data 
points. Each face is characterized by two vectors: (1) the edge vector, which connects the data points of the cells 
sharing the face and (2) the directed-area vector, which is normal to the face and has the amplitude equal to the face 
area. For each cell/face combination, the vectors are directed outward. Local grid skewness is measured at a face as 
the skew angle between the edge vector and the directed-area vector; small skew angles correspond to low skewness. 
On grids with convex cells, the maximum skewness is bounded by 90°; on more general grids with non-convex cells, 
the skew angle can approach 180°. As shown further in the paper, the grid skewness plays a critical role in defining 
accuracy and stability of FVD solutions. 

For most of the computational tests, the random node perturbation in each dimension is defined as \rh, where 
r € [—1, 1] is a random number and h is the local meshsize along the given dimension. With these perturbations, 
triangular cells in the rectangular geometry are allowed to collapse, i.e., a cell may become a zero-volume cell, albeit 
with a probability approaching zero. The random perturbations are introduced independently on all grids in grid 
refinement implying that grids of types (IV)-(VI) are grids with discontinuous metrics, e.g., ratios of neighboring cell 
volumes and face areas are random on all grids and do not approach unity in the limit of grid refinement. 

III. Finite- volume discretization schemes 

The considered model problem is the Poisson equation 

A U = f, (1) 

subject to Dirichlet boundary conditions; function / is a forcing function independent of the solution. The 2D primal 
meshes generated for this study are composed of triangular and quadrilateral cells. The FVD schemes are derived 
from the integral form of a conservation law 
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where VC/ is the solution gradient, 11 is a control volume with boundary F, and n is the outward unit normal vector. The 
general FVD approach requires partitioning the domain into a set of non-overlapping control volumes and numerically 
implementing equation (2) over each control volume. 



Figure 2. Control-volume partitioning for fi nite-volume discretizations. Numbers 0 — 12 and letters A — L denote grid nodes and primal 
cell centers, respectively. The control volume for a node-centered discretization around the grid node 0 is shaded. The control volume for 
a cell-centered discretization around the cell center A is hashed. 

Cell-centered and node-centered discretizations are considered. Cell-centered discretizations assume solutions 
are defined at the centers of the primal grid cells with the primal cells serving as the control volumes. The cell center 
coordinates are typically defined as the averages of the coordinates of the cell’s vertexes. Node-centered discretizations 
assume solutions are defined at the primal mesh nodes. For node-centered schemes, control volumes are constructed 
around the mesh nodes by the median-dual partition: the centers of primal cells are connected with the midpoints of 
the surrounding faces. These non-overlapping control volumes cover the entire computational domain and compose a 
mesh that is dual to the primal mesh. Both cell-centered and node-centered control-volume partitions are illustrated in 
Figure 2. 

A. Cell-centered FVD schemes 

In cell-centered discretizations, the conservation law (2) is enforced on control volumes that are primary cells. The 
flux at a face is computed as the inner product of the solution gradient at the face and the directed-area vector. With 
reference to Figure 2, the gradient, V//m, at the face linking nodes 0 and 4 is computed to satisfy two relations: 

1. (Vt/o 4 • e^ d J e ) = Ff ~ U N , where e e AB e = i r f ~ r ^ i is the edge \AB\ unit vector; U a and Ub are the solutions 
defined at the cell centers; r c A and r B are the cell-center coordinate vectors. 

2. (V/7o 4 • eQ 4 ° e ) = |VC/^ oce |, where = . r * r Q, is the face unit vector and |V[/^ ace | is the solution 

| r 4 -r 0 | 

derivative computed along the face [0, 4]; r" is a coordinate vector of the node i. 

The cell-centered FVD schemes considered in this paper differ only in computing | \7(jf acf: | . 

1. Cell-centered node-averaging FVD schemes 

In the node-averaging (NA) schemes, the solution derivative along the face is computed as the divided difference 
between the solution values reconstructed at the nodes from the surrounding cell centers. With respect to Figure 2, the 
solution at the node 0 is reconstructed by averaging solutions defined at the cell centers A, B. and C. The solution 
reconstruction proposed in 1318 and used in 6 is an averaging procedure that is based on a constrained optimization 
to satisfy some Laplacian properties. The scheme is second-order accurate and stable when the coefficients of the 
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introduced pseudo-Laplacian operator are close to 1 . It has been shown in 10 that this averaging procedure is equivalent 
to an unweighted least-square linear fit. The gradient is resolved from the derivative along the face and the derivative 
along the edge connecting the cell centers across the face. For the face [0,4], the face derivative is computed using 
solutions reconstructed by node averaging at the nodes 0 and 4; the edge derivative is computed using solutions at cell 
centers A and B. 

On highly stretched and deformed grids, some coefficients of the pseudo-Laplacian may become negative or larger 
than 2, which has a detrimental effect on stability and robustness. 2,3 Holmes and Connell 13 proposed to enforce 
stability by clipping the coefficients between 0 and 2. The NA schemes with clipping represent a current standard in 
practical computational fluid dynamics (CFD) for applications involving cell-centered finite volume formulations. As 
shown further in the paper, clipping seriously degrades the solution accuracy. 

The computational stencils of NA schemes are very large, so direct relaxation of the full linearization operator 
is prohibitively expensive. Defect-correction iterations are widely used as a practical solution method. A common 
driver used for defect-correction iterations is an edge-based thin-shear-layer (TSL) FVD scheme. Although the TSL 
schemes are not accurate on general grids (accuracy of TSL FVD schemes is zeroth order on general grids), the 
defect-correction method is efficient on grids with low skewing. As shown in Section IX, efficiency and stability of 
defect-correction methods with TSL drivers degrade for grids with high skewing. 

2. Cell-centered face-least-square scheme 

An alternative cell-centered scheme relies on a face-based least-square method. First, a gradient is reconstructed within 
a face using a least-square procedure. The derivative along the face is computed by inner product of the reconstructed 
face gradient with the face tangent (unit) vector. 

In this paper, we consider two approaches to determine stencils for the least-square linear fit. The stencil of type 
(A) includes the two prime cells sharing the face and their face neighbors that share at least one of the face nodes. A 
more compact stencil of type (B) typically involves the two prime cells and two auxiliary cells; one for each prime 
cell. The auxiliary cell is chosen from the pool of the cells sharing the nodes of the face as the cell closest to the prime 
cell, but not its face neighbor. The stencil of type (B) is important for discretizations on high-aspect-ratio grids of 
types (I)-(III) to represent correctly the direction of the strong coupling. 

Three types of least-square methods are considered. In Cartesian coordinates, both weighted (FWLSQ) and un- 
weighted (FULSQ) least-square methods are considered. In the FWLSQ method, the contributions to the minimized 
functional are weighted with weights inversely proportional to the distance from the face center; in the FULSQ method, 
all contributions are equally weighted. For the gradient approximation in regions with curvature, an additional approx- 
imate mapping (FAMLSQ) method is introduced. The method uses the distance function that is defined as the distance 
to the nearest boundary and normally available in practical schemes. 

B. Node-centered FVD scheme 



Figure 3. Illustration of gradient reconstruction for viscous terms on mixed grids with median-dual partition. 


The second-order accurate node-centered FVD scheme illustrated by Figure 3 represents a standard CFD approach 
to node-centered viscous discretizations. The scheme approximates the integral flux through the dual faces adjacent to 
the edge [Pq, Pf\ as 


VP • hdT = VU R ■ n R + • n L . 


( 3 ) 


ABC 
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The gradient is reconstructed separately at each dual face as follows. 

For the triangular element contribution, the gradient is determined from a Green-Gauss evaluation at the primal- 
grid element, 

VC/ i =W 0 i4. (4) 


The gradient overbar denotes a gradient evaluated by the Green-Gauss formula on the primal cell identified by the 
point subscripts. With fully-triangular elements, the formulation is equivalent to a Galerkin finite element scheme 
with a linear basis function. 1 - 2 On triangular grids of types (II) and (III) in rectangular geometries, the formulation 
recovers the five-point Laplacian stencil of the type (I) grids, independent of aspect ratio. 

For the quadrilateral element contribution, the gradient is evaluated as 

\7U r = W0234 + [|Vf/ B | - W0234 • e 04 ] e 04 , (5) 


where 


|V£H = 


U A -U 0 


l r 4 - r 5i 

is the edge gradient computed at the median Pb of the edge \I] h P\], Ui and r 6 7 ' are the solution and the coordinate 
vector of the node T), and 

r p - r p 

_ _ r 4 r 0 / 7 \ 

e ° 4 — I_p -PI G) 

l r 4 r ol 

is the unit vector aligned with the edge [ Pq , P 4 ]. Note that for grids with dual faces perpendicular to the edges, the 
edge-gradient, \VU E \, is the only contributor. This approach to the gradient reconstruction is used to decrease the 
scheme susceptibility to odd-even decoupling. 9 - 11 It has been shown in 4 - 21 that the scheme possesses second-order 
accuracy for viscous fluxes on general mixed-element grids. 


(6) 


IV. Complexity of discretization stencils 

The size of the stencil for the viscous discretization is examined for 2D and 3D cell-centered and node-centered 
FVD schemes. Estimates are made for Cartesian meshes split into triangular and tetrahedral elements and neglecting 
any boundary effects. Estimates are compared to a numerical calculation on an actual 3D turbulent viscous grid that 
include boundary effects. 

In two dimensions, two splittings of the Cartesian grid are considered. The first splits each quadrilateral cell with 
a diagonal oriented in the same direction. The second splits the cells with diagonals of face-adjacent quadrilaterals 
oriented in the opposite direction. The second splitting is slightly more analogous to the 3D splitting. In three dimen- 
sions, half of the grid nodes have 18 incident edges (32 incident tetrahedra) and half have 6 incident edges (8 incident 
tetrahedra). Each of the tetrahedra interior to an originally-hexahedral cell is defined by four nodes, each with 18 
incident edges. Each of the four surrounding tetrahedra within an originally-hexahedral cell is defined by three nodes 
with 18 incident edges and 1 node with 6 incident edges. For reference. Table 1 shows the average and maximum 
number of edges, n e d ge , connecting to a grid node. The average number of connecting edges sets the stencil size for 
the node-centered scheme as s = n e d ge + 1. The number of connecting edges is also an important factor in the NA 
schemes because the values in the control volumes faces are computed from cell-centered data averaged to the grid 
nodes. 7 


Dimension 

Hedge (Average) 

n ed ge (Maximum) 

2D 

6 

8 

3D 

12 

18 


Table 1. Edges connecting to a grid node in the split Cartesian grids. 

Table 2 shows stencil-size estimates for 2D triangular grids. There is a slight difference in the estimates from the 
two splittings (entries separated by slashes in the table), depending on the diagonalization pattern. The NA stencil is 
the largest. The face least-square (LSQ) stencil is estimated for the approach (A). 

Table 3 shows stencil-size estimates for 3D tetrahedral grids. The NA stencil is again largest. The face least-square 
stencil is only slightly larger than the stencil of the node -based discretization, in both estimation and computation. The 
discretization complexity of face LSQ stencils of type (B) is even smaller. 
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Node-centered 

Cell-centered 

NA 

Cell-centered 
face LSQ 

7 

13/16 

10/9 


Table 2. Size of the viscous stencil for 2D discretizations with triangular elements. 



Node-centered 

Cell-centered 

NA 

Cell-centered 
face LSQ 

Estimate 

13 

79 

15 

Numerical 

14 

69 

15 


Table 3. Average size of the viscous stencil for 3D discretizations with tetrahedral elements. 


V. Analysis methods 


A. Method of manufactured solution 

Accuracy of FVD schemes is analyzed for known exact or manufactured solutions. The forcing function and boundary 
values are found by substituting this solution into the Poisson equation with Dirichlet boundary conditions. The 
discrete forcing function is defined at the data points. 

1. Discretization error 

The main accuracy measure is the discretization error , E d , which is defined as the difference between the exact discrete 
solution, U h , of the discretized equations (2) and the exact continuous solution, U, to the differential equation (1) 

E d = U- U h ; (8) 


U is sampled at data points. 


2. Truncation error 

Another accuracy measure commonly used in computations is truncation error. Truncation error, E t , characterizes 
the accuracy of approximating the differential equation (1). For finite differences, it is defined as the residual obtained 
after substituting the exact solution U into the discretized differential equations. 12 For FVD schemes, the traditional 
truncation error is usually defined from the time-dependent standpoint. 20,22 In the steady-state limit, it is defined fe.g., 
in 7 8 ) as the residual computed after substituting U into the normalized discrete equations (2), 




(f (VC/ • n) dT 


J 


where |fi| is the measure of the control volume. 


(9) 


|fi| = JJ dii, (10) 

n 

f h is an approximation of the forcing function / on Q, and the integrals are computed according to some quadrature 
formulas. Note that convergence of truncation errors is expected to show the order property only on regular grids; on 
irregular grids, it has been long known that the design-order discretization-error convergence can be achieved even 
when truncation errors exhibit a lower-order convergence or, in some cases, do not converge at all. 8, 14, 15 
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3. Accuracy of gradient reconstruction 


Yet another important accuracy measure is the accuracy of gradient approximation at a control-volume face. For 
second-order convergence of discretization errors, the gradient is usually required to be approximated with at least 
first order. For each face, accuracy of the gradient is evaluated by comparing the reconstructed gradient, V r , with the 
exact gradient, V eX act, computed at the face center. The accuracy of gradient reconstruction is measured as the relative 
gradient error: 



( 11 ) 


where functions e and G define face-wise amplitudes of the gradient error and the exact gradient, respectively, 


e = \y r U h - V exact U 1 1 and G = |V eX actt/|; 


( 12 ) 


U and U h are a differentiable manufactured solution and its discrete representation (usually injection) on a given grid, 
respectively; j| • j| is a norm of interest computed over the entire computational domain. 


B. Downscaling test 

Recent publications 4 21 introduced a novel practical computational analysis method, downscaling (DS) test, for ana- 
lyzing convergence of discretization errors on irregular grids. The DS test consists of a series of inexpensive compu- 
tational experiments that account for local properties of the studied scheme; it is designed to provide estimates for the 
convergence orders of the discretization errors by comparing errors obtained on different scales. Analysis methods 
predicting the performance of DS tests have also been developed. While rigorous proofs of estimates for convergence 
of FVD truncation and discretization errors in realistic computations are still out of reach, the DS test provides an 
accurate and efficient accuracy assessment on various irregular grids, for various types of solutions and boundary 
conditions. 


C. Computational analysis of discretization errors 

On irregular grids derived from regular underlying grids by introducing local irregularities, the total discretization 
error observed in grid refinement computations is composed of (1) the smooth discretization error and (2) the local 
discretization error. The smooth discretization error is the error that would be observed in grid refinement on the 
background regular grids. An estimate of the convergence order of the smooth discretization errors can be obtained 
from observing convergence of averaged truncation errors on the background regular grids. The local discretization 
error is caused by grid irregularities. Its convergence can be directly analyzed in DS tests. For a manufactured 
solution, the discretization error profile observed on a computational grid is an accurate indicator of which component 
is dominant: large smooth error component indicates that the main error source is the solution approximation on the 
background regular grid; a spiky error profile is an indication that a major part of the discretization error is introduced 
by local grid irregularities. 


VI. Isotropic irregular grids with discontinuous metrics 

A. Grids and manufactured solution 

A sequence of consistently refined grids of either types (IV) or (VI) is generated on the unit square [0, 1] x [0, 1], 
Irregularities are introduced at each grid independently, so the grid metrics remain discontinuous on all the grids. The 
ratio of areas of neighboring faces can be as large as because a control volume can be arbitrary small, the ratio of 
the neighboring volumes can be arbitrary high. Isotropic grids randomly generated for this study have 0.01% of cell 
volumes smaller than ^-2^; the total number of grids nodes is ( N + 1) 2 ; about 0.006% of interior faces have the skew 
angle larger than 70°. The effective meshsize, h e , is computed as the L \ norm of square roots of the control volumes. 

B. Gradient reconstruction accuracy 

The accuracy of gradient reconstruction for isotropic irregular grids is first order for all methods, 5 which is suffi- 
cient for second-order discretization accuracy. As an example, the gradient reconstruction tests are performed for the 
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Figure 4. Accuracy of gradient reconstruction on isotropic irregular grids with discontinuous metrics 


manufactured solution U = sin ( 772 ; + 2 Try). For both node-centered and cell-centered formulations, gradients are 
reconstructed within the control-volume face. Figure 4 shows convergence of the L ^ norm of relative gradient errors 
on a sequence of refined grids of type (IV). Only cell-centered schemes are shown. The gradients are reconstructed 
by the FULSQ and FWLSQ methods and by the NA method without clipping. Face-based least-square methods use 
stencils of type (A). 

C. Convergence of truncation and discretization error 

The numerical tests evaluating convergence of truncation and discretization errors are performed with Dirichlet bound- 
ary conditions specified from the manufactured solution U = sin (ttx + 2 Try) and imposed on all cells linked to the 
boundary. Figure 5 shows convergence of the norms of truncation and discretization errors for the three cell- 
centered formulations on grids of type (IV). As predicted in, 4,21 truncation errors do not converge on irregular grids in 
any norm. Discretization errors converge with second order for all the three formulations considered. 




(a) Truncation Errors 


(b) Discretization Errors 


Figure 5. Convergence of truncation and discretization errors on type (IV) random triangular grids with discontinuous metrics. 

Another series of tests is performed to compare the relative accuracy of cell-centered and node-centered schemes 
on the same consistently refined mixed-element grids of type (VI). For mixed-element grids, approximately half 
of the cells are randomly split. Dirichlet boundary conditions are specified from the manufactured solution U = 
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Figure 6. Discretization error convergence on type (VI) mixed-element grids with discontinuous metrics. Prefi x CC denotes cell-centered 
scheme; NC denotes the node-centered scheme. CC-NA scheme without clipping is used. 


sin (2nx + 2i ty). Figure 6 demonstrates convergence of the Li-norm of the discretization errors. Note that on these 
mixed-element grids the number of the cells is about 50% larger than the number of nodes. The discretization errors 
of the cell-centered and node-centered FVD schemes are almost over-plotted, indicating a similar accuracy per degree 
of freedom, which implies better accuracy for the cell-centered FVD schemes on the same grids, as expected. 

D. Effects of clipping 

The tests reported in this section are performed for NA schemes and demonstrate detrimental effects of clipping on 
accuracy of gradient approximation and on the discretization accuracy. The accuracy is evaluated for the manufactured 
solution U = sin(27rt/). Dirichlet boundary conditions over-specified from the manufactured solution are imposed at 
all control volumes connecting to the boundary. Considered irregular grids of type (IV) with discontinuous metrics 
are derived from underlying isotropic (unit aspect ratio) Cartesian grids covering the unit square. Figure 7(a) shows 
an example of an isotropic random triangular grid of type (IV) with 17 2 nodes; nodes where clipping occurs in 
reconstruction are circled; about 7% of the interior nodes are clipped. 

1. Gradient approximation accuracy 

The computational gradients are evaluated within interior faces. At each interior face, the normal and tangential com- 
ponents of the computational gradient are computed and compared with the exact-gradient components at the face 
center. Figure 7(b) shows the maximum norm of the deviation between the computational and the exact-gradient 
components. The results indicate that the face gradient computed by the NA scheme with clipping does not approx- 
imate the exact gradient. The solution of the NA scheme without clipping provides a first-order accurate gradient 
approximation. Although not shown, the FULSQ and FWLSQ schemes also provide first-order accurate face-gradient 
approximations. 

2. Discretization errors 

Figure 8 exhibits convergence of discretization errors. The no-clipping NA scheme demonstrates second-order con- 
vergence; convergence of the NA scheme with clipping appears as second order on the coarse grids, but then degrades 
to the zeroth order. Although not shown, the norms of discretization errors converge with the same orders as the 
corresponding L t norms. Inspection of the discretization error profile shows that clipping greatly affects the smooth- 
error component. Plots of the discretization-error surfaces sliced along the y direction at x i — are shown in 
Figures 9(a) and 9(b) for a grid with 257 2 nodes; the figures display discretization errors at the centers of the cells 
crossed by the slice versus the y coordinate of the cell centers. For the scheme with clipping, the smooth error is 
dominant. In fact, for the NA scheme with clipping, the results of a DS test with the focal point at (x,y) = (0.3, 0.3) 
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Face gradient approximation 




(a) Random triangular grid with 17 2 nodes. Clipped nodes are (b) Relative gradient error (Loo norm) 

circled. 


Figure 7. Gradient approximation accuracy of NA schemes on isotropic irregular triangular grids with discontinuous metrics 



i 


10 10 
effective meshsize, h 


Figure 8. Convergence of NA schemes on isotropic irregular triangular grids with discontinuous metrics 


(Figure 10) show that the local errors converge at least with the first order. For the NA scheme without clipping, 
amplitudes of the smooth and local error components are comparable; both converge with the second order. 

VII. Anisotropic irregular grids with discontinuous metrics 

In this section, we study FVD schemes on irregular stretched grids generated on rectangular domains. Specifically, 
a sequence of consistently refined stretched grids is generated on the unit square [0, 1] x [0, 1] in the following 3 steps. 

1. A background regular rectangular grid with (iV+1) 2 nodes and the horizontal mesh spacing h x = jt is stretched 
toward the horizontal line y = 0.5. The coordinates of the horizontal grid lines in the top half of the domain are 
defined as 

V$+i= 0-5; y j =y j - 1 +h y p j -^ +1 ), j = y + 2, j + 3 , . . . , N, N + 1. (13) 

Here h y = ^ is the minimal mesh spacing between the vertical lines; A = 1000 is fixed maximal aspect ratio; 
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Discretization error 




(a) Node averaging scheme with clipping (b) Node averaging scheme without clipping 

Figure 9. Discretization error on a grid with 257 2 nodes along the vertical line x 



Figure 10. Downscaling test: discretization errors for NA schemes 


(3 is a stretching factor, which is found from the condition y s+\ = 1. The stretching in the bottom half of the 
domain is defined analogously. 

2. Irregularities are introduced by random shifts of nodes in the vertical and horizontal directions. The vertical 
shift is defined as A yj = ^pmin(/i^ -1 , h 3 y ), where p is a random number between —1 and 1, and h ^ _1 and 
h J y are vertical mesh spacings on the background stretched mesh around the grid node. The horizontal shift is 
introduced analogously, A Xi = j^ph x . With this random node perturbations, all perturbed quadrilateral cells 
are convex. 

3. Each perturbed quadrilateral is randomly triangulated with one of the two diagonal choices; each choice has 
probability of one half. 

An irregular random stretched grid with 65 2 nodes is presented in Figure 11(a). The manufactured solution is U = 
sin (ttx + 2?ry + . 

A recent study 5 assessed accuracy of gradient approximation on various irregular grids with high aspect ratio 
A = -jp 2> 1 . Summary of the relevant results for grids of types (I)-(VI) is presented in Table 4. The study indicates 
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10' 1 

effective meshsize, h 


10 ° 


(a) Random triangular stretched grid with 65 2 nodes 


(b) Discretization error ( L \ norm) 


Figure 11. Convergence of cell-centered FVD schemes on stretched irregular triangular grids with discontinuous metrics 


that for high-aspect-ratio grids and manufactured solutions significantly varying in the direction of large mesh spacing, 
the face-gradient reconstruction may produce extremely large relative errors (the relative errors are 0(Ah x )). 

Table 4. Relative error of gradient reconstruction 


Grids 

(I) 

(II) 

(HI) 

(IV) -(VI) 

NC 

O(hl) 

o(hl) 

0{h x ) 

0(Ah x ) 

CC-NA 

o(hl) 

O(Ahl) 

0(Ah x ) 

0(Ah x ) 

CC-FULSQ 

o{hl) 

O(Ahl) 

0(Ah x ) 

0(Ah x ) 

CC-FWLSQ 

O(K) 

0(Ah 2 x ) 

0(Ah x ) 

0(Ah x ) 


A poor gradient reconstmction accuracy, however, does not necessarily imply large discretization error. Mavriplis 16 
reported (second-order) accurate node-centered solutions even on grids with large gradient reconstruction errors. Here, 
we observe similar results for cell-centered formulations. Convergence of the L \ norms of discretization errors on a 
sequence of consistently refined stretched grids with maximum aspect ratio A = 1000 is shown in Figure 11(b). 
Although somewhat erratic, all discretization errors are relatively small and converge in grid refinement with second 
order on average. 


VIII. Grids with curvature and high aspect ratio 

In this section, we discuss accuracy of FVD schemes on grids with large mesh deformations induced by a combi- 
nation of curvature and high aspect ratio. The grid nodes are generated from a cylindrical mapping where (r, 6) denote 
polar coordinates with spacings of h r and hg, respectively; the innermost radius is r = R. The grid aspect ratio is 
defined as the ratio of mesh sizes in the circumferential and the radial directions, A = The mesh deformation is 
characterized by the parameter F: 


r _ -R(l - cos (hg)) _ Rhg _ A h e 
h r 2h r 2 

The following assumptions are made about the range of parameters: R = 1, A 1, and \'h r <C 1, which implies 
that both h r and hg are small. For a given value of A, the parameter F may vary: r > 1 corresponds to meshes 
with large curvature-induced deformation; f <C 1 indicates meshes that are locally (almost) Cartesian. In a mesh 
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refinement that keeps A fixed, F = O(Ahg) asymptotes to zero. This property implies that on fine enough grids with 
fixed curvature and aspect ratio the discretization error convergence is expected to be the same as on similar grids 
generated on rectangular domains with no curvature. 

We are focusing on difficulties specific to high-T grids with large curvature-induced deformations. The difficulties 
arise for functions that are predominantly varying in the radial direction of small mesh spacing. For functions predom- 
inantly varying in the direction of small mesh spacing, gradient reconstruction in rectangular geometries is expected 
to be accurate. 5 

Three types of 2D grids are considered for the cylindrical geometry: (I) regular quadrilateral grids ; (II) regular 
triangular grids derived from the regular quadrilateral grids by the same diagonal splitting of each quadrangle; (III) 
random triangular grids, in which regular quadrangles are split by randomly chosen diagonals, each diagonal orienta- 
tion occurring with probability of half. Grids of types (I) and (II) are regular; grids of type (III) are irregular because 
there is no periodic connectivity pattern. The grid types are shown in Figures 12. For a better visualization of the 
grid topology, low-F cylindrical grids are shown. Random node perturbation is not applied to high-T cylindrical grids 
because even small perturbations in the circumferential direction may lead to non-physical control volumes. 




(a) Type (I): Regular quadrilateral grid (b) Type (II): Regular triangular grid 


(c) Type (III): Random triangular grid 


Figure 12. Typical cylindrical grids. 


A. Accuracy of gradient approximation 

It was observed and confirmed analytically, 16- 17- 19 that the unweighted-least-square gradient approximation is zeroth 
order accurate on deformed grids with high V. Least square method with inverse-distance weighting has been proposed 
and shown to provide an accurate gradient approximation at nodes for node-centered formulations. 

The study 5 shows that both FULSQ and FWLSQ methods suffer accuracy degradation on high-T grids. To im- 
prove the accuracy of gradient approximation, a least-square minimization in a mapped domain has been proposed. 
Two mapped least-square methods have been considered: an exact mapping (FEMLSQ) method that uses the polar co- 
ordinates directly; and a more general approximate mapping (FAMLSQ) method that is based on the distance function, 
defined as the distance to the nearest boundary, normally available in practical schemes. Applicability of FEMLSQ is 
limited to model problems with analytical boundary shape. 

The more general FAMLSQ method approximates the FEMLSQ method by applying the least-square minimization 
in a locally constructed coordinate system. 

/’■ = /o + «£ / + A»/. (15) 

The local coordinates, (£', ?/), are constructed using the distance function, which provides information on the closest 
boundary point. First the distance function is reconstructed at the face center. The coordinate vectors at the face 
center are defined as a unit //'-directional vector pointing in the direction opposite to the closest boundary point and 
its orthonormal ^'-directional vector. The //'-coordinate is its distance from the boundary and the ^'-coordinate is the 
projection of the vector connecting the stencil point with the face center onto the ^'-direction. For the cylindrical 
geometry, the distance function is the shifted radial function, r — R. The gradient approximation accuracy for high-F 
grids of types (I)-(III) has been studied in 5 and is summarized in Table 5 . The least-square gradient reconstruction 
methods shown in the table use stencils of type (A), but the results are the same for stencils of type (B). 
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Table 5. High-T grids: relative errors of F-gradient reconstruction 



(I) 

(II) 

(III) 

NC 

o{hl) 

0(h g ) 

0(hg) 

CC-FULSQ 

0{hl) 

0(1) 

0(1) 

CC-FWLSQ 

o(h 2 e ) 

0(1) 

0(1) 

CC-FEMLSQ 

o(h 2 r ) 

o(h 2 ) 

o(h 2 ) 

CC-FAMLSQ 

o(h 2 ) 

o(h 2 ) 

0(hl) 

CC-NA 

o(h 2 e ) 

o(h 2 e ) 

O(Ahg) 


B. Discretization error convergence 

Discretization errors of cell-centered schemes are compared with the errors of the node-centered scheme on a sequence 
of stretched high-r grids of type (III). The tests are performed for the manufactured solution U = sin(7rr). The 
computational grids are derived from background regular cylindrical grids with the radial extent of 1 < r < 1.04 
and the angular extent of 10°. The background grids have twice as many nodes in the radial direction than in the 
circumferential direction. The stretching performed in the radial direction provides a fixed maximal aspect ratio 
A = 21817 throughout the grid refinement study. The maximal value of parameter T changes approximately from 240 
to 30. The stretching ratio is changing as (3 = 1.94, 1.38, 1.175, and 1.083. A representative 17 x 33 grid is shown in 
Figure 13(a). 



— e- 

- CC-FULSQ 

B — 

- CC-FWLSQ 

0— 

- CC-FAMLSQ 

— v — 

NC 



Effective Mesh Size 


(a) High-T grid 


(b) Discretization errors 


Figure 13. Comparison of cell-centered and node-centered discretizations on high-I grids of type (III) 

Convergence of the L \ - norm of the discretization errors is shown in Figure 13(b). The cell-centered (CC) least- 
square gradient reconstruction methods use stencils of type (A). The errors of the node-centered (NC) solution and 
the CC-FAMLSQ solution converge with second order and are almost over-plotted on fine grids indicating the same 
accuracy per degree of freedom. The errors of other two cell-centered least-square methods are significantly higher 
and appear convergent with first order. Inspection of the discretization stencils reveals that the NC stencil correctly 
reflects the strong coupling in the direction of small mesh spacing; the CC-FULSQ and CC-FWLSQ stencils fail to 
identify the strong coupling direction. 

Least-square methods with compact stencils of type (B) for gradient reconstruction are better suited for high- 
aspect-ratio grids of types (II) and (III). Figure 14(b) shows convergence of cell-centered methods on a sequence of 
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consistently refined stretched grids of type (III) with the radial extent 1 < r < 1.001 and the angular extent of 10°. The 
fixed maximum aspect ratio is A = 100, 000; F « 550 on the coarsest 9 x 17 grid and T « 70 on the finest 65 x 129. 
The stretching ratio is changing as (3 = 1.57, 1.25, 1.11, and 1.055. The 17 x 33 grid is shown in Figure 14(a). The 
manufactured solution U = sin(27r r + 7t/6). 


1.001 r 

1 .0005 - 
1 - 

0.9995 - 
0.999 - 
0.9985 - 
0.998 - 
0.9975 - 
0.997 - 
0.9965 - 

0.996 1 1 1 1 1 1 1 1 1 1 

-0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1 




r=R(1-cos(h 0 ))/h r 


(a) High-r grid 


(b) Discretization errors 


Figure 14. Cell-centered discretizations on high-r grid of type (III) 

Four cell-centered schemes are tested: the FULSQ scheme with compact stencils of type (B), the FAMLSQ 
schemes with stencils of types (A) and (B), and the NA scheme. Similar to the results of the previous tests (Fig- 
ure 13(b)), the discretization errors of the FAMLSQ scheme with stencils of type (A) are the smallest. However, the 
errors of the FULSQ scheme with stencils of type (B) are not much larger, showing significant improvement over the 
least-square schemes with stencils of type (A). The discretization errors of the FAMLSQ scheme with stencils of type 
(B ) are over-plotted with the FULSQ errors, indicating that the FULSQ convergence is optimal for the given stencil 
choice. The discretization errors of the NA scheme are several orders of magnitude higher than the errors of other 
schemes. 


IX. Stability of defect-correction iterations for node-averaging schemes 

In this section, we apply the local mode Fourier (LMF) analysis to analyze stability of defect-correction iterations 
on regular skewed grids of type (II). A sketch of a skewed grid is shown in Figure 15(a). The analysis is performed 
for a discrete Fourier test function e^ exlx+er > lr i\ where (0 X , 0 V ) are normalized Fourier frequencies in the horizon- 
tal x-direction and in the skewed ^-direction, respectively. Figure 15(a) illustrates the notation. For node-centered 
formulations, the Fourier symbol of the discrete operators is a scalar functional of the normalized frequencies; for 
cell-centered formulations, the symbol is a 2 x 2 matrix that couples Fourier components defined at the low and upper 
triangles. 

Defect-correction iterations with the TSL driver are applied to the target NA scheme without clipping on grids with 
high skewing. The grid skew angle is defined as the maximum skew angle over all faces. The stability of iterations 
enhanced by the main-diagonal contributions from pseudo-time stepping characterized by the viscous CFL number. 
The results of the LMF analysis shown in Figure 15(b) demonstrate the relation between the maximum grid skew angle 
and the CFL number required for stability of the defect-correction iterations. The CFL number drops to two on grids 
with high skewing, effectively making the scheme explicit. The predictions of the LMF analysis have been precisely 
confirmed by numerical testing. Note that the CFL requirements for stability of defect-correction iterations for the NA 
scheme with clipping are exactly the same. 
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(a) Skewed grid of type (II) 


(b) CFL vs. skew angle 


Figure 15. Maximum viscous CFL number providing stability of defect-correction iterations on skewed grids of type (II). The driver 
operator is edge-based TSL scheme, the target operator is NA scheme without clipping. 


X. Conclusions 

Node-centered and cell-centered schemes are compared for finite-volume discretization of Poisson’s equation as a 
model of the viscous flow terms. The comparisons are made for two classes of tests: the first class is representative of 
adaptive-grid simulations and involves irregular grids with discontinuous metrics; the second class is representative of 
high-Reynolds number turbulent flow simulations over a curved body. 

Among the compared schemes, grids, and solutions, the node-centered scheme looks the most attractive. The 
scheme is known to deteriorate to first order for inviscid solutions on irregular mixed-element grids , 4,21 but for viscous 
terms, its accuracy is always second order and comparable with the accuracy of the best cell-centered schemes at 
equivalent number of degrees of freedom; also the linearization stencil of the node-centered scheme involves the least 
number of neighbors. 

The cell-centered node-averaging (NA) scheme has three significant drawbacks. The first is that the NA scheme 
with clipping can fail to converge to the exact solution. Note, however, that the clipping is introduced mainly for 
stability of the inviscid solution and can be eliminated for the viscous terms. The second drawback is that the scheme 
demonstrates inferior accuracy on certain triangulated highly anisotropic grids in curved geometries. The third draw- 
back is that the scheme’s linearization stencil is quite large in comparison to the other schemes. Only the edge-based 
linearization terms are typically used for iterative solutions in a defect-correction setting. For high grid skewing, un- 
avoidable in highly anisotropic triangular/tetrahedral grids used for high Reynolds number simulations, the scheme 
has explicit stability limitations. 

The cell-centered schemes using least-square face gradient reconstruction have more compact stencils with a com- 
plexity similar to the complexity of the node-centered scheme. For the first class of tests representative of adaptive-grid 
simulations, the face least-square methods, the NA method without clipping, and the node-centered method demon- 
strate second-order convergence of discretization errors with very similar accuracies per degree of freedom. Accurate 
solutions have been obtained even for combinations of grid/manufactured-solution/scheme for which gradients are 
approximated with large relative errors. For simulations on highly anisotropic curved grids, the least-square methods 
have to be amended either by introducing a local mapping of the surface anisotropy or modifying the scheme stencil 
to reflect the direction of strong coupling. 
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